Virtually all dynamic systems can be modeled so that signal processing techniques can be applied to the models to extract useful information about the system. Many diverse signal processing techniques have been developed with an eye towards ever increasing the speed of signal extraction, and increasing the ability to deal with many different kinds of signals, some of which when modeled are not well behaved.
Signal processing systems and methods find wide application in a plethora of fields and industries. Magnetic resonance imaging (MRI), radar and sonar imaging, medical imaging, surface science and electronics, non-destructive testing, sensor applications, ultrasound, acoustical reception, geophysics, and other fields all use some form of signal processing to analyze signals gathered from a signal source which are then further conditioned to be operated upon by a signal processing system having appropriate means to process the data. It will be recognized by those with skill in the art that signal processing routines are usually accomplished on standard digital computers which are interfaced to a data gathering system or data input system that provides the signal information to the computer for analysis.
Signal processing methods which have heretofore been developed can be generally broken down into the following discrete steps: gathering signal data from a source or object under examination, converting the signal data to a form usable by a digital computer, inputting the converted data to a digital computer for analysis, and operating on the converted data with a particular data processing method which characterizes the data to the desired appropriate form. The data from the source may be input in many different ways such as, for example, an energy spectrum, a dispersion relationship, an intensity spectrum, or an amplitude versus frequency plot. Whatever the form of the data input to the computer, the basic step of operating on the data according to a particular signal processing method is always performed to produce meaningful output information about the object or source.
There are probably as many signal processing techniques as there are signals to be analyzed. Those skilled in the art will be readily familiar with the fast Fourier transform, the LaPlace transform, and many different numerical techniques which qualify as signal processing methods. A very useful class of signal processing methods requires that the signal received be first approximated to an n-th order linear difference equation. Signal processing systems which utilize such an approximation require that signal information be modeled as a system function which is a ratio of polynomials and which can be completely determined by the pole-zero distribution of the system function in the complex plane. When the signal of interest gathered from the signal source or object to be identified involves noise signals, the problem of detecting possibly damped, sinusoidal signals in the noise involves very high degree polynomials which are usually derived by fitting a linear recursion sequence to the signal samples or to estimated autocorrelation samples.
In this situation, the polynomial's degree is typically one-quarter to one-half the number of sample points, and so polynomial orders of several hundred are common. The polynomial coefficients are generally small, indicating well-conditioned polynomials, that is small perturbations in the coefficients cause similar small changes in the root locations in the complex plane. In the past, polynomial coefficients have been derived from noisy data, and small measurement errors by the data gathering devices do not change the underlying system response in a drastic manner. In these types of systems, the roots of the polynomials are usually evenly distributed in angles close to the unit circle which is a plot of the roots in the complex plane, although it is possible that a few of the roots may substantially deviate from the modulus one of the circle. Sinusoidal components of the signal are derived by identifying the signal roots, and spurious non-signal roots are considered as an autoregressive approximation to the noise in the signal.
There are many ways to characterize a signal as an n-th order linear difference equation. One of the most common techniques for approximating a signal as a linear difference equation is to apply Prony's method to the signal. See, e.g., S. Lawrence Marple, Jr., Digital Spectral Analysis with Applications, Chapter 11, pp. 303-349 (1987), the teachings of which are specifically incorporated herein by reference. Prony's methods fits a deterministic exponential model to the data to obtain a spectral interpretation of the energy spectral density. According to Prony's method when as many data samples are used as there are exponential parameters, an exact exponential fit to the data may be made. This relationship can be expressed as follows: ##EQU1## There are thus p equations in Equation 1 wherein 1&lt;n&lt;p, which can be expressed in matrix format as: ##EQU2## This matrix represents a set of linear simultaneous equations that can be solved for an unknown vector of complex amplitudes.
According to Equation 1, it is apparent that a solution to some homogeneous linear constant coefficient difference equation may be found. In order to find this difference equation, it is desired to define a polynomial .phi.(z) that has z.sub.k exponents as its roots, which may be expressed as: ##EQU3## Equation 3 can be expanded into a power series so that the polynomial can be represented in a summation such as: ##EQU4## Equation 4 has complex coefficients a[m] such that a[o] equals one. The polynomial in Equation 4 is generally denoted as a "characteristic equation" associated with the linear difference equation.
Prony's method requires fitting p exponentials to 2p data samples and can be summarized in three steps. The first step is to obtain the polynomial coefficients associated with Equation 4. The second step is to calculate the roots of the polynomial which are defined by Equation 4. Many types of polynomial factoring routines have been devised to perform this second step such as, for example, the CPOLY computer program taught in the above-referenced Marple text. In calculating the roots of the polynomial, a damping coefficient, .alpha..sub.i, and a sinusoidal frequency f.sub.i, are determined from the roots z.sub.i of Equation 4 after index shifting and coefficient substitutions for convenience using the relationships: EQU .alpha..sub.i =ln.vertline.z.sub.i .vertline./T sec.sup.-1 ( 5) EQU f.sub.i =tan.sup.-1 [Im{z.sub.i }]/Re{z.sub.i }2.pi.T Hz (6)
To complete Prony's procedure, the roots so computed are used to construct the matrix of elements which can then be solved for the p complex parameters h[1], . . . h[p] of Equation 2. The amplitude A.sub.i and the initial phase .theta..sub.i may be determined from each h.sub.i parameter, according to the relationships: EQU A.sub.i =.vertline.h.sub.i.vertline. ( 7)
and, EQU .theta..sub.i =tan.sup.-1 [Im{h.sub.i }/Re{h.sub.i }] radians (8)
The foregoing explanation of the generalized Prony's method assumes that the number of data points, N, does not exceed the minimum number needed to fit a model of p exponentials, that is, N is greater than 2p. However, in the real world the number of data points usually exceeds the minimum number, and the system is then said to be "overdetermined." In such a case, it has been necessary to approximate numerically the polynomial roots. Numerical factorization schemes of a very high degree of well-conditioned polynomials have been developed to accomplish signal processing in this situation. See, e.g., K. Steiglitz and B. Dickinson, "Phase Unwrapping by Factorization", IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. Ass. P-30, No. 6 (December 1982), the teaching of which are specifically incorporated herein by reference. Other numerical schemes have been developed with an eye towards solving massive root polynomials.
While these numerical techniques have proven themselves useful for certain signal processing problems, numerical determination of polynomial roots is a difficult problem at high values of order n. Search techniques based on Newton-Raphson iterations often encounter numerical overflow or underflow of the high order polynomials and their derivatives everywhere except in the vicinity of the unit circle, thereby creating errors in the root solutions. Root solving techniques generally resort to "deflation," which is the process of repeatedly dividing the polynomial by the roots that these techniques can locate. Deflation tends to produce generally ill-conditioned intermediate polynomials having a huge dynamic range in their coefficients. The deflated polynomials can exhibit chaotic behavior, wherein small errors in the coefficients cause massive shifts in the root positions, thereby also producing spurious root results.
In the worst case scenario, errors in deflated polynomials that are accumulated during deflation will leave the final polynomials, which would normally contain the actual signal roots, without useful precision. Root polishing techniques which have been employed in conjunction with deflation methods may recover some of the precision, but are generally defeated by the inability to evaluate the original polynomial in a sufficiently large convergence area around a root which is far from the unit circle.
Other numerical techniques have been developed and employed to approximate and model systems in light of these problems. For example, least squares analyses have been used with Prony's method of approximating a signal spectrum, and the aforementioned phase unwrapping by factorization according to the Steiglitz et al. article has also been developed. These numerical techniques, and others, are available for approximating solutions to the linear difference equation constants.
However, a basic problem inherent in all the numerical methods mentioned above is encountered when determining the coefficients of a polynomial with a digital computer having appropriate software, that is, the problem of overflow and underflow of computer registers which invariably occurs during processing of the computer software associated with the methods. Consider the simple case of a polynomial of order 200, P(z)=z.sup.200 -1. Using Digital Equipment Company's VAX d-float format having a 55-bit mantissa and an 8-bit exponent, the CPOLY root solver taught in the Marple text returns serious errors in 15 percent of the root positions on the unit circle. Furthermore, using the search-deflate-polish technique with a worst case intermediate-order polynomial of N=100 where roots lie in one complex half-plane, with a 106-place mantissa representation of the roots approximately one-half of the mantissa is lost in d-float format. Extending precision defers this problem to higher orders without resolving the overflow and underflow problems. Standard software implementations of CPOLY thus usually limit the acceptable polynomial order to n=50, for example. This problem is inherent in virtually every prior numerical-type signal processing method or system.
Therefore, prior numerical techniques for calculating the coefficients of a polynomial to perform signal processing simply fail with polynomials of order greater than about 50. While the above-referenced techniques of signal processing perform nominally for limited applications, it is desirable to develop new methods and systems of signal processing which provide highly accurate root solutions of characteristic polynomials, provide accurate calculation of polynomial coefficients, and are reliable for high orders. Furthermore, the root solving techniques which are used during object-identification signal processing methods should be easily implementable on standard, digital computers and adaptable for many types of signals and signal processing applications. The prior numerical techniques referenced above simply do not satisfy these criteria.